* Estimates levels of the DV for 2 groups, gp1 and gp2
program define levels2, rclass
   version 16.1
	syntax varlist [if] [pweight], GP1(varname) GP2(varname) ///
	   [LABEL1(str) LABEL2(str) LABEL3(str) OPTS(str) ALLFIVE NODRAW]
	marksample touse
	
	if "`weight'" ~= "" {
	   local wt [`weight'`exp']
   }		
	if "`label1'"~="" & "`label2'"~="" {  // labels for legend, if desired
	   local label1 label(`label1')
		local label2 label(`label2')		
	}
	else {
	   local legend legend(off)
   }
	
	forvalues i=1/2 {
	   table axn5 if `touse' & `gp`i''==1 `wt', c(mean `varlist' n `varlist') format(%4.2f)
	}
	
	tempname lvl1 lvl2
	local gopts scheme(s1mono) plotregion(style(none)) mlabel format(%3.0f) `legend'
	
	* Distinguish all five action-cost combinations
	if "`allfive'" ~= "" {
	   forvalues i=1/2 {
		   reg `varlist' i.axn5 if `touse' & `gp`i''==1 `wt', robust
		   margins i.axn5
		   mat `lvl`i'' = r(table)
      }				
	   if "`nodraw'" == "" {
			coefplot ///
				(matrix(`lvl1'[1,]), ci("`lvl1'[5,] `lvl1'[6,]") msymbol(O) mlabpos(2) `label1')  ///
				(matrix(`lvl2'[1,]), ci("`lvl2'[5,] `lvl2'[6,]") msymbol(Oh) mlabpos(10) `label2'), ///
				`gopts' `opts'
      }				
   }			
	else {
	   forvalues i = 1/2 {
		
			* Estimates if no action
			reg `varlist' if `touse' & axn5==0 & `gp`i''==1 `wt', robust // mean if no action
			margins
			mat `lvl`i'' = r(table)
			
			* 5% reduction in emissions. Average over prices
			reg `varlist' i.axn5 if `touse' & (axn5==1 | axn5==2) & `gp`i''==1 `wt', robust
			margins, asbalanced // equal weight over prices
			mat `lvl`i'' = `lvl`i'', r(table)
			
			* 25% reduction in emissions. Average over prices.
			reg `varlist' i.axn5 if `touse' & (axn5==3 | axn5==4) & `gp`i''==1 `wt', robust
			margins, asbalanced // equal weight over prices
			mat `lvl`i'' = `lvl`i'', r(table)
						
			mat colnames `lvl`i'' = none low high						
			* mat list `lvl`i''
		}
		if "`nodraw'" == "" {
			coefplot ///
				(matrix(`lvl1'[1,]), ci("`lvl1'[5,] `lvl1'[6,]") msymbol(O) mlabpos(2) `label1')  ///
				(matrix(`lvl2'[1,]), ci("`lvl2'[5,] `lvl2'[6,]") msymbol(Oh) mlabpos(10) `label2'), ///
				coeflabel(none = "No Action" low = "Cut 5%" high = "Cut 25%") ///
				`gopts' `opts'			
         }				
	}
	
	di _n(2) "Levels if `gp1'"
	mat list `lvl1', noheader format(%5.2f)
	di _n(2) "Levels if `gp2'"
	mat list `lvl2', noheader format(%5.2f)
	return matrix lvl1 = `lvl1'
	return matrix lvl2 = `lvl2'
		
end
